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INTRODUCTION 


The method described below provides an easy and 
inexpensive way to evaluate the Fourier series coeffi- 
cients for an empirical function. The computation is 
carried out by taking the values of the function f(x) at 
certain values of x multiplied alternately by +1 or -1 
and followed by summation. No weighting coefficients 
are needed, which facilitates the computation. The 
points x, at which the value of the function is needed, 
however, vary from one harmonic to another. There- 
fore it is preferable to plot the function at a certain 
scale and use a transparent template for finding these 
points for each harmonic. 

In all numerical methods for Fourier analysis, the 
higher harmonics have a disturbing effect on the results. 
This is displayed here so that the coefficient of the n:th 
harmonic contains all the harmonics 6 nm +1. (m = 1, 2, 
3...). If, however, the whole Fourier spectrum is calcu- 
lated, the results are easily corrected for these higher 
harmonics, this in particular if the series converge 
rapidly. 


THEORY 


The real and imaginary parts of the exponential Fou- 
rier transform F (h) of function f (x) 


w 
-@ 
(1) 
| +2 ihx 
f dh = [F(-hy] 
can be expressed as 
Re F (h) = Fe (h) = cos 2 dx 
-@ 


(2) 


@ 
Im F = Fo(h) =- f(x) sin 2 hx dx 


where f. and .. are the even and odd parts of f. 
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(3) 


| | fe - =f (-x) 


The two following operators are also needed in the.deri- 
vation 1) 


oo 
rep, f (x) = > f (x - » t) 
v=-@ 


(4) 


oo 
comb, f (x) = » 3 f(v t) 6 (x- » t) =f (x) rep, 6 (x) 


v=-@ 


6 (x) is the Dirac deltafunction and t is a constant para- 
meter. Operators (4), when transformed, yield two 
important relationships 1) 


f (x)] comb, F (h) 
T 


(5) 


*[comb, f (x)] = repy F (h) 
t 


Same equations hold also for the inverse transformation. 
In (5) let t = 1/n and f(x) = f,(x) or f(x) = f(x). We obtain 
when both sides are divided through by n 


RY [ rep; fe (x)] = comb, Fe (h) 
n 


oo 
> Fe (ay) (h-ny) 
vy =-@ 


(6) 


rep; fy (x)] =i comb, Fo (h) = iF, (nv) 6 (h-n 
n 


Obviously spectra (6) are discrete with 46-peaks at each 

multiple of n. If expressions similar to (6), but t = 1/2 n, 

are written down, they have dé-peaks only at even ’ 
multiples of n. The difference of these two expressions 

thus leaves the odd multiples of n. 


= [2 repy fe (x) f, (x)] = > (nv) 6 (h-nv) 
n n v odd 


(7) 


= [2 rep, fy (x) “rep (x)] = iF, (nv) 6 (h-ny) 
n 2n v odd 


Using the transformation 


(8) = 6 (h-nv) 


and properties 


(9) F, (h) = Fe (-h); Fo (h) = (-h) 


for a real function f(x) equations (7) give 


‘ 
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Be (x)= [2 rep, fe (x) “rep fe (x) ] =2 > (nv) cos 2 % nvx 
n 


2n 
odd 


(10) ; 


Bo (x)= [2 rePy fy (x) 1 fo (x)] = Fy (nv) sin2 nvx 
odd 


Let further p(x) be a periodic function with period equal 
to unity and let us denote by f(x) the function that is equal 
to p(x) over one period and zero elsewhere. 


when -$< x<} 
(11) 


when x <-4 or x>$4 


Now the Fourier series coefficients a, and b, of p(x) are 


i 

a, =2 cos2 wnxdx =2 [fe cos 2 wn x ax 
-4 
4 


b, = 2 posin 2% nx dx =2 sin 2 nx dx 
-4 -@ 


(12) 


where f. and f. are now the even and odd parts of (11). 
By equations (2) we have 


a, = 2F, (n) 
(13) n>o 


9 
’ 
| 
| by = ~ 2 Fy (n) 
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and finally from (10) and (13) 


| Be (xX) = a, cos 2 MN X + ag, Cos2 3nx + a5, cos2 SNX+... 


(14) 


Bo (x) = bp sin 2% nx +bg, sin 2 3nx sin 2 W5nx+... 


Finding suitable values for x in g, and g, we can still 
eliminate all coefficients whose subscripts are divisible 
by 3n, These turn out to be x = 1/12n for g, and x= 1/6n 
for g,- Thus 


(15) 


1 3 


By equations (10) and (4) equations (15) can be worked out 
to yield 


@ . 
1 
v=-@ 


(16) 


co 


The range of summation from -o@ to oo is only formal, 
for by definition (11) f(x) is zero outside the interval : 
(-1/2, 1/2) and the sums contain actyally only 2n terms. 

According to equations (16) the coefficients a, and | 
by of successive harmonics of periodic function p(x) are 
calculated in the following way: We take one period of the 
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function according to (11) and find the even and odd parts 
of it (3). Next we find the points (1-6 »)/12n or(1-3r)6n 
respectively and evaluate the sums (16) for n= 1, 2, 3»-N 
until the series converge to a value small enough in 
comparison with the accuracy required. And finally the 
coefficients obtained in this way can be corrected for the 
higher harmonics 5n, 7n, l1ln, etc. 

Even the zero order harmonic a, can be evaluated by 
this method. This can be accomplished in two different 
ways: Obviously the value of p(x) or f(x) at the origin is 


(17) p(o) =f(0) =4a,+ ay 


If the cosine series foefficients a, (n>0) have been 
evaluated in the manner suggested above, we can easily 
compute the sum 


Hence 


(18) a, (0) - ay] 
v=1 


Another way of finding a, is suggested by the first of 
equations (6). 
By an inversion 


‘ 

> ay- 

v=1 
@ 

| 


Fe (nv) cos 2a nvx 
=1 


+ fe (x) = Fe (0) + 2 
n 


2 


(19) 


=ta,+ any Cos 
v=1 


If, within the limits of required accuracy, a, = 0 for 
| 


@ 
1 
= 


or 


@ 
_ 2 


=-@ 


Keeping in mind equations (3) and (11), according to which 
f, (-1/2) = f, (1/2) = 4p (1/2), expression (20) turns out 
to be the trapezoidal rule for the integral of f, the spac- : 
ings of ordinates being equal to 1/n. 
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TEMPLATE 


The main drawback of the method is that the points 
(1-6»)/12n or (1-3»)/6n, at which the values of the 
function are needed, do not generally coincide. This 
difficulty can be overcome by using a template engraved 
on a transparent sheet (lucite, celluloid, etc). 

Two sheets, each having the breadth of one period, 
which has been standardized to, say, 40 cm, are divided 
by straight lines into horizontal strips. One of the sheets 
is used to the calculation of a,, the other for bi On the 
first strip (corresponding to n= 1) of the a, sheet are 
marked the points (1-6 )/12 with vertical lines extending 
across the strip, and the same procedure is repeated for 
each n on both sheets. The sign of each term is indicated 
by colouring the lines corresponding to positive (» even) 
red and the negative (» odd) black. Now, when f, (x) and 
f, (x) have been plotted on a plotting paper, the ordinates 
can be read easily by moving the templates up and down 
on the paper. 

By using properties f, (x)=f, (-x) and f, (x) = - f, (-x) 
it suffices to make templates covering only half of the 
period. The points from the negative half of the period 
are shifted on the positive half with proper sign. On co- 
sine template the signs are not altered, on sine template 
they are opposite. Fig. 1 depicts cosine and sine tem- 
plates for n=1-7. 
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Fig. 1. Sine and cosine template for coefficients a, and b,, n=1 - 7. 


Terms of positive sign have been drawn with continuous lines and 


those of negative sign with dashed lines. 
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INVERSION 


In many instances, e.g. crystal analysis, we know 
the Fourier coefficients (a,, b,) and want to obtain the 
function p (x) or one period of it f (x). This problem can 
also be handled by the template method. 

By equations (13) we have the real and imaginary part 
of the Fourier transform of f (x) at integral values n of h 


F, (h) =R, F (hy =4a n =h = integer > 0 
(21) e e h 


Fo (h) = Im F (h) = - $y (by = 0) 


And, because ReF (h) is an even and ImF (h) an odd 
function, functions F, and F, can be extended to the nega- 
tive values of h, too. 


Fe (-h) = 4 ap, n =h = integer >o 
(22) | 


Fy (-h) =4 by 


It is easily verified that the periodic function p (x) is 
an infinite repetition of one of its periods 
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(23) p(x) = repy f(x) 


but there is an infinite number of functions f, (x)* which, 
when repeated at intervals 1, give the same function. 


(24) p(x) = rep, f, (x) 


Among these f(x) is the "simplest" in the sense that 
the repetitions never operlap one another. Obviously both 
functions (23) and (24) must have the same series, or by 
(21) and (22) the transforms of f and f, coincide at h = 
integer. At nonintegral values they may differ from one 
another. It follows from this property that by joining the 
points F (n) by an arbitrary (continuous) curve we get the 
transform Fy, (h) of one of functions f, (x). 

One of the basic properties of Fourier transforms is 


(25) = [f, =F, @mh) m = constant 


or by the second of equations (1) 


(26) [Fy (-2mh)] = fy 


Function F), (-2 mh) is aperiodic but becomes periodic : 
by operation | 


(27) P(h) = rep, (-2mh) = P, (h) (h) 


| 
| 
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In analogy with (24) the transform of P (h) coincides 
with that of F(-2mh) at points x = integer. Or by (13) 


fk (ga) =2a, 


fk, -2by = integer > 
(28) 
-v 

‘he (Gm) =? 


fk (se) =2b, 


where a, andb, are the cos and sin series coefficients 
of P, (h) and P, (h) and f,,, and f,,. are the even and odd 
parts of f,. Finally 


a,-b, 
(29) (gm) = 2 

a,+by v¥<o 
and by (24) 
(30) = > fy = fy 


In this way 2m points of p (x) are obtained to the 
period. It can be shown by a detailed discussion that if 
the coefficients a,, b,, are known up to the N:th harmonic, 
the best accuracy in" p (x) is obtained by taking m= N. 
A higher value of m gives more points to the period but 
they do not improve the detail in p (x). If on the other 


4 
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hand m < N, the points »/2m are too far apart giving a 
poorer detail and some of the information is lost. 

The degree of overlapping of functions f, in (30) de- 
pends essentially on the smoothness of the curve by which 
the points F(2nm) are joined. The smoother curve gives 
less overlap. Theoretically, if the interpolation is done 
with Hadamard’s 2) interpolation formula 


sin 2 7 m (h-n) 
2am (h-n) 


(31) F,(2mh) = F(2nm) 


no overlap occurs 3). In practice a smooth hand-drawn 
curve approximates very well this function. 

In practice, then, the inversion is done by dividing 
the half period into N equal parts (m = N) and plotting the 
coefficients a, and b, at the points of division. Then the 
points are respectively joined with a smooth curve and 
analysed with the template. The possible overlapping is 
corrected by equation (30). 

Theoretically there are no limits as to how high har- 
monics can be analyzed in this way, but the work needed 
for each harmonic grows in proportion to n, since there 
are 2n terms to add. 

A detailed analysis indicates that if all the ordinates 
in (16) are read with a certain accuracy + A, the error 
in a, and b, is approximately + A/Vn. This means that 
the higher harmonics are more accurate than the lower 
ones in the absolute sense. 
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